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Abstract 

The  application  of  uncertainty  analysis  to  curve  and  surface  fits 
obtained  by  the  method  of  Least  Squares  is  described.  The 
primary  obstacle  to  the  routine  implementation  of  the  uncertainty 
analysis  procedure:  the  derivation  of  the  sensitivity  derivatives,  has 
been  removed.  Analytic  expressions  for  the  derivatives  applicable 
to  fits  of  arbitrary  order  have  been  derived,  and  a  step-by-step 
procedure  for  their  incorporation  within  a  computer  program 
provided.  A  review  of  the  techniques  for  the  construction  of  curve 
and  surface  fits  of  arbitrary  order  is  included.  The  use  of 
uncertainty  analysis  as  an  aid  to  the  assessment  of  the  suitability  of 
a  particular  fit  is  demonstrated  by  the  application  of  these 
procedures  to  a  series  of  examples  employing  a  generic  data  set 
typical  of  experimentally  derived  data.  Some  results  include: 
quantification  of  the  penalty  for  using  a  higher-order  fit  when  it  is 
not  appropriate,  and  the  fact  that  reduction  of  uncertainty  in  the 
data  to  be  fitted  is  more  effective  at  reducing  the  uncertainty  in  fit 
coefficients  and  fitted  values  than  simply  increasing  the  amount  of 
data  used  to  construct  the  fit.  A  computer  program  implementing 
these  procedures  is  available  from  the  authors. 

Administrative  Information 

This  work  was  sponsored  by  the  Office  of  Naval  Research  (Code  333)  under  Contract 
N0001497WX30267  and  Program  Element  601 153N. 

Introduction 

The  use  of  the  method  of  Least  Squares  for  the  construction  of  curve  fits  as  models  for  experimental 
behavior  is,  for  better  or  worse,  a  universal  practice.  At  an  early  stage  of  training  in  engineering 
disciplines  and  in  the  sciences,  one  is  introduced  to  the  procedure  and  quickly  learns  to  produce  linear 
and  quadratic  fits.  Less  well  known  is  the  general  implementation  of  the  procedure  for  polynomial  fits 
of  arbitrary  order  and  the  ability  to  replace  the  usual  set  of  basis  functions  with  alternatives  which  may 
be  more  appropriate.  A  useful  example  of  the  latter  approach  is  the  use  of  products  of  polynomials  of 
arbitrary  order  for  the  construction  of  a  surface  fit  to  a  function  of  two  independent  variables.  The 
indiscriminate  or  incorrect  use  of  the  method  can  lead  to  faulty  conclusions  when  modeling  behavior, 
and  because  the  technique  is  so  popular  and  straightforward  to  apply,  the  potential  for  misuse  is  large. 

Measures  to  assess  the  goodness  of  the  fit  exist  and  serve  as  a  check  on  the  suitability  of  a  particular 
model,  but  these  measures  can  be  insensitive  and,  considered  in  isolation,  misleading.  One  must  also 
consider  the  degree  to  which  measurement  uncertainty  existing  in  the  data  to  be  fitted  will  propagate  into 
the  fit  coefficients  and  in  fitted  values  obtained  from  the  use  of  the  model  when  judging  the 
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appropriateness  of  a  particular  fit.  This  paper  discusses  the  application  of  uncertainty  analysis  to  curve 
and  surface  fits  of  arbitrary  order  determined  using  the  method  of  Least  Squares.  The  general  procedure 
requires  the  calculation  of  various  partial  derivatives,  termed  sensitivity  derivatives,  which  can  be  quite 
formidable  if  attempted  in  a  brute-force  fashion.  Instead,  analytic  expressions  for  the  computation  of 
these  derivatives  when  the  method  is  applied  to  curve  and  surface  fits  have  been  found;  these  formulas 
may  be  implemented  within  a  computer  program,  thereby  making  the  use  of  uncertainty  analysis  an 
integral  part  of  the  procedure  for  the  construction  and  assessment  of  an  arbitrary  order  fit.  The  manner 
in  which  these  derivatives  are  determined  is  linked  to  the  procedure  for  calculating  arbitrary  order  fits; 
for  this  reason,  the  general  method  for  constructing  curve  and  surface  fits  will  be  reviewed. 

The  calculation  of  uncertainty  associated  with  coefficients  for  and  with  fitted  values  from  a  Least 
Squares  fit  is  attributed  to  Coleman  and  Steele  and  may  be  found  in  the  new  edition  of  their  text  ‘.  As  of 
this  writing,  the  text  is  scheduled  to  be  released  in  February  1999.  Because  the  method  does  not  appear 
in  their  previous  text  and  has  been  available  only  through  private  courses  offered  in  the  past  few  years, 
the  uncertainty  analysis  procedure  along  with  general  instructions  for  its  use  have  been  reproduced  here 
for  convenience.  Direct  calculation  of  the  sensitivity  derivatives  for  the  linear  and  quadratic  curve  fit 
cases  is  provided  to  allow  an  independent  check  on  the  implementation  of  the  general  method.  A  step- 
by-step  outline  for  incorporating  the  general  procedure  within  a  computer  program  is  included.  All  of 
the  techniques  described  here  have  been  executed  in  a  FORTRAN  computer  program  which  is  available 
upon  request  from  the  authors. 

Finally,  the  use  of  uncertainty  analysis  as  an  aid  to  the  assessment  of  the  suitability  of  a  particular  fit  is 
demonstrated  by  a  series  of  simple  examples  employing  a  generic  data  set  which  is  typical  of 
experimentally  derived  data.  The  results  illustrate  the  general  use  of  the  method,  quantify  certain 
previously  held  suspicions  and  introduce  some  surprises.  The  next  section  begins  with  a  description  of 
the  general  procedure  for  the  derivation  of  an  arbitrary  order  fit. 


Computation  of  an  Arbitrary  Order  Least  Squares  Curve  Fit 

Given  n  pairs  of  experimental  data  (x;  ,y, )  for  i  =  1, 2, . . . ,  n ,  one  can  fit  an  mth  order  polynomial  to  the 
data  of  the  form: 


y  —  Cj  +  CjX  +  + . . .  +  Cm+j  x  1  ■t  j 

k= 0 

where  the  choice  of  m  is  constrained  such  that  n>m  + 1 .  The  values  of  the  unknown  coefficients  are 
determined  by  the  criterion  used  to  choose  the  best  fit  to  the  data.  The  Least  Squares  criterion  chooses 

the  coefficients  of  the  fit  such  that  the  sum  of  the  squares  of  the  distances  of  the  measured  data  (x,  ,y,) 

from  the  fitted  data  (x,. ,y,)  is  a  minimum.  In  other  words,  we  seek  to  minimize  the  quantity 
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Because  the  (x,  ,>>,)  are  known,  q  will  depend  only  upon  the  coefficients  ck  for  k  =  1, 2, + 1.  To 
minimize  q,  we  compute  partial  derivatives  with  respect  to  each  coefficient  and  set  each  resulting 
equation  equal  to  zero  to  obtain 
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Rearranging  yields  the  following  set  of  m  + 1  simultaneous  equations: 
c,n  +  c2  Z  x;  +  c3  Z  x(  + . . .  +  cm+)  Zx<  -  Zt/ 

c,  X + c2  Z  * 2 + c3  Z  + •  •  • + c*+i  Z  <+1 = Z 
c,Zx,2  +c2Zt3 + c3  Z  x/4 + ••• +c«+iZxr2 = Zx.2t 

c,Zxr  +c2Zxr’  +  c3  Z  *r2 + •  •  • + ^,+1  Z  *2m = Z  *r  t, 

These  equations  can  be  easily  represented  in  matrix  form  and  are  known  as  the  normal  equations. 
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Here  we  have  a  (w+ 1)  x  (m+ 1)  matrix  of  known  quantities  multiplying  an  unknown  (m  +  l)  x  1 
coefficient  vector  to  obtain  a  known  (m  + 1)  x  1  right-hand-side  vector.  There  are  a  variety  of  methods 
available  for  the  solution  of  a  simultaneous  system  of  linear  equations.  Note  that  the  {m  + 1)  x  (m  + 1) 
matrix  above  is  symmetric  and  positive-definite  which  implies  that  very  efficient  and  fast  simultaneous 
equation  solvers  may  be  employed  for  the  solution.  Alternatively,  if  the  matrix  is  nearly  singular,  then 
singular  value  decomposition  techniques  may  be  applied.  These  procedures  are  described  in  Numerical 
Recipes3. 

To  build  the  above  set  of  equations  within  a  computer  code,  one  begins  by  forming  the  n  x  (m+ 1) 
matrix  C  and  its  (m  + 1)  x  n  transpose  CT  : 
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Then,  the  normal  equations  defined  in  Eqs.  5  can  be  easily  developed  using  matrix  multiplication  as 
follows: 


CTC  c  =  CTy  ,  (7) 

and  the  known  matrices  CTC  and  CTy  can  be  passed  to  a  simultaneous  equations  solver  to  solve  for  the 
unknown  coefficient  vector  c. 

Following  Holman4,  quantities  which  describe  the  degree  of  goodness  of  the  fit  are  the  standard  error  of 
the  fitted  values  given  by 
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and  the  correlation  coefficient,  R,  which  is  computed  from 


R  = 


where  the  standard  deviation  ofy, 


1 

n  - 1 


t(y, -yf 


(9) 


Note  that  some  texts  omit  the -  coefficients  in  the  numerator  and  denominator  of  the  expression 

n— ... 

that  defines  R.  The  coefficient  of  determination,  used  by  some  texts,  is  also  used  as  a  measure  of  the 
goodness  of  the  fit  and  is  simply  the  square  of  the  correlation  coefficient.  This  quantity  is  usually 

defined  such  that  the  — - —  coefficients  in  the  numerator  and  denominator  are  omitted. 
n  — ... 


Uncertainty  Method 

Given  n  pairs  of  experimental  data  (x,.,y,.),  along  with  associated  bias  errors,  (Bx).  and  (^).,  and 

precision  errors,  (Px).  and  (pJ  ,  the  method  due  to  Coleman  and  Steele 1  determines  how  these  errors 

propagate  through  the  governing  equations  to  produce  uncertainty  in  the  resulting  fit  coefficients  and  in 
future  values  predicted  from  the  fit.  Consider  the  simple  case  of  a  linear  fit 
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The  equations  for  the  coefficients  are  the  solutions  to  Eqs.  5  for  the  case  m  =  1  and  are  given  by 


i=\  i=\  /= 1  /=! 


^Lx,y, 


M  ^  /=! 


and  c-,  = 
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To  compute  uncertainty  in  fit  coefficients,  the  method  treats  equations  such  as  Eqs.  11  (or  Eqs.  21 
below)  as  propagation  equations  and  then  applies  standard  techniques  for  the  propagation  of  uncertainty. 
These  equations  are  presented  in  the  next  section.  Then,  combining  Eqs.  10  and  11  one  can  track  the 
uncertainty  propagated  into  future  fitted  values  which  will  be  described  in  the  following  section. 


Uncertainty  in  Curve  Fit  Coefficients 

Denote  the  bias  uncertainty  propagated  into  the  kfh  coefficient  ck  as  (Bc  ) ;  this  quantity  may  be 
computed  from 


where  we  have  assumed  correlated  biases  among  the  xj ,  among  the  yi  and  also  between  x,  and  y, .  If 
any  of  these  biases  among  the  raw  data  values  can  be  shown  to  be  uncorrelated,  then  the  corresponding 
terms  in  Eq.  12  should  be  excluded. 

The  precision  uncertainty  propagated  into  the  coefficient  ck ,  (Pc)  ,  is  obtained  from  a  similar  equation; 

however,  since  precision  errors  are  considered  to  be  random,  all  precision  errors  among  the  raw  data 
values  are  uncorrelated. 


To  arrive  at  an  overall  uncertainty,  ( Uc)  ,  for  the  coefficient  ck ,  one  must  combine  the  bias  uncertainty 
computed  in  Eq.  12  with  the  precision  uncertainty  obtained  in  Eq.  13.  There  are  two  generally  accepted 
methods  for  computing  (uj)  .  The  first  method  is  denoted  as  the  root-sum-square  (RSS)  method  and 
the  total  uncertainty  is  given  by 


This  uncertainty  is  considered  to  be  a  95%  coverage  estimate  when  B  and  P  are  95%  confidence  values. 
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The  second  method  is  to  combine  bias  and  precision  uncertainties  by  simply  adding  them  (ADD 
method): 

dXdB\dpX-  (15> 

This  method  produces  a  total  uncertainty  which  provides  about  99%  coverage  when  B  and  P  are  95% 
confidence  values  and  when  neither  B  nor  P  is  negligible  relative  to  the  other.  However,  when  either  B 
or  P  is  negligible  relative  to  the  other,  clearly  the  total  uncertainty  cannot  be  better  than  a  95% 
confidence  estimate.  To  carry  out  these  computations,  we  need  expressions  for  the  partial  derivatives 
(sensitivity  derivatives)  that  appear  in  Eqs.  12  and  13.  For  the  simple  case  of  a  linear  fit,  the  derivatives 
are  applied  to  Eqs.  11.  We  make  use  of  the  formula  for  the  derivative  of  a  quotient,  and  for  the  linear 
case,  we  have: 
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For  the  case  of  a  quadratic  fit,  the  equations  for  the  coefficients  are  the  solutions  to  Eqs.  5  for  the  case 
m  =  2  and  are  given  by: 
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;z< 

j+Z*f[Z*Z  *?-«Z*?]+Z*.‘ 

»z*,2-(z*,)2 

Z4Z*,‘Z*-Z*,2Zv, 

FZ*?. 

yi[Z*iZ*i2-»Z*i 

TZ*. 

f[«Zx/T 

,-Zx/Zt,-. 

Z^jZ^Z^-G*2) 

Z^2[Zx/Z^/  -  Z*2  Z  y] 

2]+Z*/[2>‘2>< -"Z^l+Z*,2 

+Z*.![Z*iZ>’.  -"EwM1: 'y 

«Z<  - 

;[«Z*.2  - 

(Z-.T 

-(Z*,2 

l 

(21) 

)2] 

z*,2[z*,z*Mz*2)2; 

+Z4 

£x,.]£x*-«£x,3]  +  2>,4  n 

Z*,2  -G 

-1  \2 

W  J 
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Applying  derivatives  to  these  expressions  is  a  tedious  task.  After  much  algebra,  we  find: 


dnumx 

dxt 

=  Tjx\3x?'Lx?y 

1 

4^ 

X 

M 

X 

+Zx-2[3A2Zx/T 

+  4.v;Zt,  -2x,.y,Z^2  - 

4^Zx2T,] 

+Z x?  [y.  Z  xf +  2x,  Z  x-y>  -  6xf  Z  y> 

+2*,j',Zxi  +Z  x?y] 

+Z^,4[2x,Z  yt  - 

Z  XM  -  Ti  Z  xi  ] 

dnun\ 

dy, 

=  Z*/[-*;Z*2] 

Z*,2- 

dnum2 
d  xt 

=Z*i[2*/Z*iV, 

,-4x,3Zt,]  +  Zx2 

2x,y,] 

£*,+3  -*£*=] 

+Zx'[2aZt,- 

2wciy\  +  Y,xi[nyl- 

Z  y] 

+Zx(T,[Zx2  ~3nxf] 

dnum2 

=Zx,2[x,2Zx-- 

x/Zx.2]+Zx>3[Zx.2-m 

fl+Z*,4  [«,-£>,] 

dnum 3 
dx. 

=  Zx.[3-T2 Zt  - 

-2*/y(Z*/]  +  Z*/2| 

[*Z*i  -4-^Zt,  +2nx,y,  +  Z*/J 

+Zx/[Zt,  -ny, 

]  +  Z  X<  T/  [2^/  Z  xi  ~ 

-3  rvcf 

|  +  Z  X/2 y>  [2nxi  _  2Z  X;  ] 

dnum2 

dy, 

=Zx,[-x,2Zx,] 

+Z*/a[*iZ*/ +m 

i-I. 

i*/2]+Z*/[Z*i-«i] 

d  den 

dx, 

=Z*.[6x2Zx.2  + 

4Xj  Z  Xi  ~  4x/3  Z  xi . 

|  +  Zx,2K-6,,Zx,2] 

+Z *1  [2Z xi  - 6nxi  ]+Zx<4 \lnx'  ~2Hx<' 


(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 


d  den 


=  0 


(29) 


These  equations  are  then  substituted  into  Eqs.  16  and  17  to  find  the  needed  derivatives.  Summarizing, 
one  computes  the  bias  and  precision  uncertainty  propagated  into  the  fit  coefficients  from  Eqs.  12  and  13 
making  use  of  Eqs.  11  and  16-20  for  the  m=  1  case  or  of  Eqs.  16-17  and  21-29  for  the  m  =  2  case.  The 
total  uncertainty  in  fit  coefficients  is  then  determined  from  Eqs.  14  or  15. 
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Uncertainty  in  Predicted  Values  from  Curve  Fit 

After  calculating  the  coefficients  of  the  fit  equation,  one  can  supply  a  new  value  X  in  order  to  compute  a 
new  value  Y using: 

Y(X)  =  c{  +c2X  +  c3X2  +  ...  +  cm+]Xm  forX  chosen  such  that  xmin  <  X<xmax  (30) 


We  wish  to  consider  the  uncertainty  that  will  propagate  into  this  fitted  Y(X).  Different  equations  are 
used  depending  upon  whether  or  not  there  is  uncertainty  associated  with  X.  The  most  general  case  is  the 
one  in  which  there  is  uncertainty  associated  with  X.  The  bias  uncertainty  that  propagates  into  the  fitted 
Y(X)  for  this  case  may  be  computed  from 


J 

\  f 


\dxJt 


W,  (B,). 


4i&)  (S)  w  w,- 4  MMi  w  to- 


(31) 


f  dY\ 

+Ux, 


tZKdX) 


d_Y_ 


b*  (*,), 


where  we  have  assumed  the  most  general  case  of  correlated  biases  among  the  x,  (and  X),  among  the  yt 
(and  X)  and  also  between  xi  and  y, .  If,  as  described  above,  any  of  these  biases  can  be  shown  to  be 
uncorrelated  then,  of  course,  the  appropriate  terms  in  Eq.  3 1  should  be  excluded. 

The  precision  uncertainty  that  propagates  into  the  fitted  Y(X)  is  obtained  from 


Now,  to  evaluate  these  expressions,  we  need  to  determine  the  partial  derivatives: 
One  applies  the  derivatives  to  Eq.  30  which,  for  the  linear  case,  may  be  written  as 


3Y  DY 
dx{  ’  dy, 


and 


3Y_ 

dX 


Y(X)  = 


n  n  n 


/= I _ /-l  M 


»=1 


(33) 
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The  needed  derivatives  are  found  to  be: 


dY  dc,  dc 7 


X  + 


dx,  dx,  +  dx, 


dY  dc,  dc2 
- L  +  — -X  + 


dyt  dy,  dy, 
dY 

-~i~X  =  ^2  T  2c3  A"  H 


dc2 

7y, 


X2 3  +...  +  ~iLXm  , 
dxi 

(34) 

S  c 

■X2  +...+— r^X"'  ,  and 

(35) 

...+mcnnXX"'-\ 

(36) 

0  c  @  c 

where  — —  and  — —  are  as  calculated  previously  for  the  linear  or  quadratic  cases  above.  Note  that 

dx-,  dy, 

3  Y  3  Yr  3  Y 

- , - and - are  functions  of  X  and  must  be  recomputed  for  each  new  value  of  X. 

dx ,  dy,  dX 


Now,  given  a  value  X,  we  can  compute  a  value  Y from  Eq.  30  and  an  associated  BY  and  PY  from  Eqs.  3 1 
and  32.  From  the  latter  two  quantities  we  can  determine  a  total  uncertainty  UY  from  Eqs.  14  or  15.  This 
process  should  be  repeated  for  Xj  which  vary  throughout  the  range  xmin  <  X}  <  xmax ,  and  a  series  of 

Y(X)j  and  (c/r)  will  result.  Then,  the  uncertainty  in  fitted  Y(X);  values  should  be  plotted  along  with 

the  fitted  curve  as  follows.  Plot  the  original  (x,  ,y, )  data  values  along  with  the  fitted  curve 

Y(X)j  =  c,  +c2Xj  +c2Xj  +...+c„„xXj  obtained  from  the  new  [Xj  ,Y^  data  pairs.  Then  plot  above 

and  below  this  fitted  curve  the  two  additional  curves:  Y{X)j  +(f'rF)/  and  Y{X) j  ~{uY)j  ■  The  latter  two 
curves  then  give  an  indication  of  the  total  uncertainty  in  fitted  Y(X)j  values  across  the  entire  range 
*min  *  Xj  <  Xmax  . 


Summary  of  Fitting  Procedure 

To  summarize,  then,  one  can  compute  the  total  uncertainty  in  the  fit  coefficients  as  well  as  in  the  fitted 
Y,  values  by  performing  the  following  steps: 

1 .  From  the  (x,  ,^(  )  data  pairs,  compute  the  fit  coefficients.  Examples  for  the  linear  and 
quadratic  cases  are  given  in  Eqs.  1 1  and  21,  respectively. 

2.  From  the  (x,  ,y,)  data  pairs,  compute  the  sensitivity  derivatives.  Examples  for  the 
linear  and  quadratic  cases  are  given  in  Eqs.  16-20  and  Eqs.  16-17,  22-29,  respectively. 

3.  Compute  uncertainty  in  fit  coefficients  (Bc)k  and  (Pe)k  from  Eqs.  12  and  13, 
respectively. 
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4.  Form  the  total  uncertainty  (Uc)k  from  Eqs.  14  or  15. 

5.  Supply  a  value  X ,  in  the  range  xmin  <  Xj  <  xmax ,  and  for  this  X}  compute  Y(X)J 

from  Eq.  30  and  4r~,  4—  and  44  from  Ecls-  34‘36- 
3xt  ay )  aX 

6.  For  this  X} ,  estimate  Bx  and  Px  . 

7.  For  this  X} ,  compute  BY  and  PY  from  Eqs.  3 1  and  32,  respectively. 

8.  For  this  Xn  form  (UY).  from  Eqs.  14  or  15. 

9.  Repeat  steps  5  to  8  for  X}  throughout  the  range  (x,  )mjn  <  XJ  <  (xi  . 

10.  Plot  the  original  (x,,y,)  data  values  along  with  the  fitted  curve 

Y(X) ,  =  Cj  +  c2Xj  +  c3Xj  +  . . .  +  cm+1XJ'  obtained  from  the  new  (x.  ,Yj)  data  pairs 

generated  in  step  9. 

11.  Then  plot  above  and  below,  respectively,  the  curve  generated  in  step  12,  the  two 
additional  curves:  Y(X)j+(uy).  and  Y{X)J-{ur)j.  The  latter  two  curves  then 
give  an  indication  of  the  total  uncertainty  in  fitted  Y(X)j  values  across  the  entire 
range  xmin  <  Xj  <  xmax . 
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Calculation  of  Sensitivity  Derivatives 

As  can  be  seen  from  the  previous  sections,  the  method  is  straightforward  to  apply;  the  primary  difficulty 
is  the  derivation  of  the  sensitivity  derivatives  for  a  fit  of  arbitrary  order,  and  we  will  show  how  this  may 
be  accomplished  in  this  section.  One  proceeds  by  first  developing  analytic  expressions  for  the 
coefficients.  This  can  be  done  by  employing  Cramer’s  rule  for  the  solution  of  the  simultaneous 
equations  in  Eqs.  5.  Cramer’s  rule  finds  the  solution  for  each  coefficient  in  c  as  a  quotient  of 
determinants;  the  denominator  is  simply  the  determinant  of  the  CTC  matrix  for  all  coefficients,  and  the 
numerator  is  the  determinant  of  the  matrix  formed  by  replacing  the  tfh  column  of  the  CTC  matrix  with 
CTy  for  coefficient  ck .  Therefore,  for  a  nonsingular  CTC  matrix  we  can  write: 


c,  = 


Z-T°T,  5>,‘ 

Zx't,  Zx,2 


x'"y. 

4  1  Si 


ilJ+1 


X*r 

X*r' 

I*,1" 


2>* 

5>; 


X*; 

X*,2 


X*; 

X*r 


1 


X«f 

Xx1  - 

Zt°t, 

X*,' 

X*?  ••• 

Zx‘t, 

X*: 

X*r‘  - 

Zxfy, 

X*,° 

X*,1 

Z*r 

X*' 

X*.2 

X*,"*1 

X*r 

X*:*1  - 

X*,2" 

(37) 


dc , 


dc , 


IX*:  X*r'  -  X*,2’ 

The  sensitivity  derivatives  that  are  needed  for  uncertainty  calculations  are  the  quantities:  and  ^ 

for  i  =  1, 2, . . . ,  n  and  k  =  1, 2, . . . ,  m + 1 .  As  an  example  of  these  calculations,  consider  the  computation 

dc. 


of 


dx* 


dc. 


X*.°  X*.v,  •••  X*; 

X*f  Xv.'v.  -  X*: 


X*r  -  X*; 


X 


X*,2" 


X,-  d  X: 


Zx? 

Zx,' 


Zx‘ 

Zx<2 


Zx' 

z< 


Zxr  Z^r1  -  Zx,2M 


(38) 


dnum  d  den 

den  — — - -num- 


dx, 


dx: 


den 2 


den 


dnum  dden 
dx ,  Ck  dx , 
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The  computation  reduces  to  the  determination  of  derivatives  of  determinants  in  the  numerator  and  the 
denominator.  To  see  how  this  may  be  accomplished  for  any  values  of  m  and  n,  we  will  first  consider  the 
simple  case  of  a  3x3  determinant. 


The  evaluation  of  a  derivative  of  the  general  case  of  an  (m+ 1)  x  (m+ 1)  determinant  can  be  carried  out 
by  computing  the  sum  of  m+l  determinants;  the  k?h  determinant  in  the  sum  is  formed  by  replacing  each 
element  in  the  k?h  column  by  its  derivative  with  all  other  columns  remaining  the  same.  This  is  an 
important  observation  because  each  of  the  determinants  in  the  sum  can  be  formed  efficiently  within  a 

dc\ 

computer  code.  For  example,  if  one  wished  to  compute  the  derivatives  — —  for  the  case  m  =  2 ,  one 

w  Xj 

would  form  the  quantities: 
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dc{  d 

Zx<° y>  Zx‘  Zx<2 
Zx‘t-  Zx»2  Zx/ 
Zx,2x  Zx/  Zx-4 

1 

x~ 

II 

Xx,°  Xx/  Zx<2 

Zx‘  Zx/2  Zx3 

Zx2  Zx3  Zx-4 

o  num 


dden 


UH.  .  -num-  .  ,  3  , 

dxi  dx,  1  dnum  dden 

den 2  den  dxi  1  dxt 

where  the  derivatives  of  the  numerator  and  denominator  would  be  found  from 

.  o  x*,1  z-2  ix-°  z xf  z xh>  z*;  2x- 

=  i*,0^  Zx,2  Zx,3  +  Z*.V/  2x'  Zx<3  +  Z*/*  Zx-2  3x2 

X‘  2x'yi  Xx.3  Zx,4  Zx.2^/  3x,2  Xx,4  Zxfr,  Zx.3  4x-3 

^  0  2>;  2>2  £x?  i-f  Zx,2  Zx,°  Zx/  2*; 

^=-1*,*  Ex?  Z*?  +  Zx,‘  2r?  Zx?+Zx‘  Ex*  Jr? 

2x?  5>?  X!x.4  Zx/2  3x,2  lx;  Xx,2  Zxf  4x<3 

.  .  .  .  dc,  „  .  .  _ 


d  den 


To  obtain  the  entire  set  of  derivatives, 


for  i  =  1, 2, . . . ,  n ,  one  must  form  the  quantities  given  in 


Eqs.  40  and  41  for  each  /  for  /  =  1, 2, . . . ,  n .  Similarly,  to  compute  -jy ,  for  the  case  m  =  2 ,  one  would 
form  the  quantities: 


Zxi°  y< 

Ex,1 

Zx<2 

ZX'T/ 

Ex,2 

Zx,3 

a  \ 

Zx,2* 

Ex? 

Zx-4 

0y, 

Zx.° 

Z x! 

Ex? 

Zx' 

Zx.2 

Ex? 

Zx.2 

Zx3 

Ex? 

dnum  dden 

den— - num— —  ,  -  -  , 

d  yi  d  yi  1  onum  dden 

den2  deny  dyt  1  dy, 

where  the  derivatives  of  the  numerator  and  denominator  would  be  found  from 


dden 
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8  num 


8  den 

7 7 


x,°  Xx'  Ix- 

Xx-V/  o  Ix,2 

Z  xi° y<  Zx‘  0 

x'  Yjxf 

+ 

0  Zx.3 

+ 

M 

J* 

vT 

M 

X 

O 

x,2  Xx»3 

Zx,2»  o  Zx< 

Z  x?y<  Zx/  0 

xi°  Zx/  Zx.2 

x!  Zx,2  Zx<3 
xf  Zx?  Z>,4 


0 

2>; 

Zx-2 

Zx-° 

0 

Z-,2 

Zxi° 

Z  x\ 

0 

0 

Zx,3 

+ 

Zx‘ 

0 

Z*,3 

4* 

Z*.1 

Z-,2 

0 

0 

s*? 

Z^4 

Zx' 

0 

Zx‘ 

Z*,! 

Z-,! 

0 

(43) 


Here  we  see  that  the  calculations  are  simpler  requiring  the  calculation  of  only  one  determinant  for  the 
numerator  (true,  for  any  order  m)  and  with  the  denominator  not  a  function  of  y, . 


Summarizing,  the  required  sensitivity  derivatives, 


<?c. 


and 


8  c,. 


for  i  =  1,2, ...  ,n  and 


8x,  “  8y, 

k  =  \,2,...,m  +  \  consist  of  two  nx(m  + 1)  arrays.  A  recipe  describing  one  possible  implementation  of 
these  calculations  within  a  computer  code  for  any  m  and  n  is  as  follows. 


1.  For  the  computations  dealing  with  coefficient  ck ,  form  a  set  of  m  +  1  arrays  of  size 
(m  + 1)  x  (m+l)  and  initially  fill  the  arrays  with  the  elements  of  the  CTC  matrix.  For 
coefficient  ck ,  replace  the  j  =  k,h  column  of  each  of  these  m  + 1  arrays  with  CTy- 

2.  Repeat  step  1  for  each  of  the  m+l  coefficients  storing  all  data;  these  arrays  will  be 

used  for  numerator  computations.  Repeat  step  1,  one  additional  time  (m  +  2'^time), 
where  each  of  the  arrays  are  filled  with  the  elements  of  the  CTC  matrix  only.  This 
latter  set  of  arrays  will  be  used  for  denominator  computations.  Note  that  all  of  this 
data  can  be  stored  within  a  four-dimensional  array  with  the  following  dimensions: 
m  +  l  rows  indexed  by  I  (where  /  is  a  different  index  than  i  which  indexes  original 
data  values),  m  +  l  columns  indexed  by  j ,  m+l  matrices  required  for  derivative 
calculations  for  each  coefficient  indexed  by  /,  and  m  +  2  sets  corresponding  to  m  +  l 
numerator  calculations  for  each  coefficient  and  one  denominator  calculation  (common 
to  all  coefficients)  indexed  by  k.  Denoting  the  four-dimensional  array  by  AIjlk  and 
looping  over  all  four  indices,  we  have 

\CTY,  j  =  k 
Aij,k  - 1 CTC,,  j*k' 

Note  that  the  k  index  can  loop  to  m  +  2  as  required  since  j  cannot  equal  k  for  this  case. 
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3.  Now,  for  the  k'h  set,  replace  the  j  =  l'h  column  in  the  Ith  matrix  in  the  set  with  — — 

C*  Xj 

of  the  previous  column  contents.  With  loops  over  7,  l  and  k  this  can  be  accomplished 
as  follows. 


0  7  =  1 


1 

J* 

1 

to 

V 

otherwise 

l  =  k 

0 

7  +  /<3 

l^k 

(7  +  /-2)  x/+/‘ 

'3  otherwise 

4.  Repeat  step  3  for  each  of  the  m  +  2  sets  indexed  by  k. 

5.  Compute  the  determinant  of  each  of  the  m+i  arrays  per  set.  Sum  these  determinants 
and  store  the  sum  in  a  one  dimensional  array  indexed  by  k. 


dnum 

6.  Repeat  step  5  for  each  of  the  m  +  2  sets.  These  sums  represent  — —  for  each  of  the 


m  + 1  coefficients  with  the  last  set  corresponding  to 


d  den 
dx. 


dx. 


.  At  this  point  all  of  the 


data  in  the  four-dimensional  array  formed  in  step  2  has  been  used,  and  we  proceed  to 

d 


the 


*y, 


calculations. 


7.  For  the  k'h  set,  replace  the  j  —  l'h  column  in  the  llh  matrix  in  the  set  of  m+ 1  arrays 
with - of  the  previous  column  contents.  One  may  use  the  four-dimensional  array 

#yt 

as  it  exists  after  step  6;  it  does  not  have  to  be  recreated.  With  loops  over  7 ,  /  and  k 
this  can  be  accomplished  as  follows. 

{0  otherwise 

*/-'  l  =  k 

8.  Repeat  step  7  for  m+ 1  sets  indexed  by  k.  Do  not  bother  with  the  m  +  2nd  set  because 
the  denominator  is  not  a  function  of  y, . 

9.  Compute  only  the  determinant  of  the  l-k'h  array  in  the  set  since  the  other 
determinants  in  the  set  will  be  zero.  Store  this  determinant  in  a  one-dimensional  array 
indexed  by  k. 
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10.  Repeat  step  9  for  m+ 1  sets  and  set  the  m  +  2  element  of  the  one-dimensional  array 


dnum 


equal  to  zero.  These  elements  represent  — 


for  each  of  the  m+  1  coefficients 


with  the  last  element  corresponding  to 


d  den 


=  0. 


dc. 


11.  Using  the  derivatives  computed  in  steps  6  and  10,  compute  the  derivatives  and 


— L  for  the  k"'  coefficient  using  Eqs.  40  and  42.  Store  the  results  in  two  two- 

dy, 

dimensional  arrays  with  dimensions  n  x  {m  + 1)  where  the  first  dimension  corresponds 
to  i  =  1,2,...,  n  original  data  values  and  the  second  to  k  =  1, 2, . . . ,  m  + 1  coefficients. 

12.  Finally,  repeat  steps  1-1 1  for  i  =  1, 2, ... , n  filling  the  arrays  described  in  step  1 1 . 

To  debug  the  implementation  one  can  independently  calculate  the  derivatives  for  the  cases  m  =  1  and 
m  =  2  by  programming  within  a  subroutine  the  derivative  expressions  derived  above  in  an  earlier 
section.  Results  obtained  from  this  routine  can  then  be  compared  to  the  results  obtained  from  the  brute- 
force  subroutine. 


Computation  of  an  Arbitrary  Order  Least  Squares  Surface  Fit 

Given  triplets  of  experimental  data  (x,  ,yj , z, )  for  i  =  1,2,...  ,n,  where  the  data  are  presumed  to  be 
related  by  a  function  of  the  form  z  =  f{x,y),  one  can  construct  a  surface  fit  with  a  functional  form 
obtained  from  the  product  of  two  order  polynomials: 


mx 

my 

z  = 

HaP^xP 

— 1 

o 

ii 

_q=0 

mx  my 

=  YYckxpy‘l , where  £  =  #  (m,.  -t-lj  +  p  +  l  , 

p=0  q= 0 


(44) 


and  there  are  (mx  +  \)(mv  + 1)  terms  in  the  sum.  Note  that  for  a  given  amount  of  data,  the  choice  of  mx 
and  my  is  constrained  by  the  condition  that  n  >  (mx  + 1  )(tny  + 1) .  When  mx  =my  =  2,  we  have: 


z  =  Ci  +c2x  +  c3x2  +c^y  +  cixy  +  cbx2y  +  c1y1  +  csxy2  +c9x2/  .  (45) 

An  advantage  to  this  formulation  is  that  one  can  choose  mx  and  my  to  be  different  values.  For  example, 
if  one  believes  the  data  to  be  linear  in  one  independent  variable  and  quadratic  in  the  other,  then  a  surface 
could  be  constructed  using  mx=  1  and  my=2.  For  even  finer  control,  an  implementation  of  this 

method  in  a  computer  program  could,  for  a  given  choice  of  mx  and  my ,  allow  the  user  to  retain  or  omit 
any  terms  in  the  desired  functional  form.  This  is  easy  to  accomplish  and  will  be  described  further  below. 

The  values  of  the  unknown  coefficients  are  determined  by  the  Least  Squares  criterion.  Therefore,  we 
seek  to  minimize  the  quantity 


16 


«  2  »  ,  v2 

tf  =  X(z/~z/)  =  X(z/  - c,  - c2x,  -  . . .  - qy,  -  cMx,y,  -  ...  - c(„„+l)(„y,+rix;">;"v )  .  (46) 

/=1  M 

Because  the  (x,  are  known,  q  will  depend  only  upon  the  coefficients  ck  for 

k  =  l,2,...,(mx  +  l)(m  + 1) .  To  minimize  q,  we  compute  partial  derivatives  with  respect  to  each 
coefficient  and  set  each  resulting  equation  equal  to  zero  to  obtain 


dq 

dci 

n 

=  2X(Z/  _C1  ~C2Xi 

1=1 

'•-W-Wi-" 

-C  Vm  VWV 

)  (—  1) = o 

dq 

dc2 

II 

to 

sMa 

1 

1 

O 

K> 

1 

/wr 

)(-x,)  =  ° 

dq 

dck 

=  2X(z/-c,-c2x,-. 
1=1 

■•-Vz-Wr- 

C(«lX+I)(w^+l)Al  ✓  / 

)(->,)  =  »  (47) 

dq 

<?ck+l 

«  , 

=  2X(z,-c,-c2x(-. 

f«l 

—  c  Xntxvmy 

C(OTJf+I)(W^+l)A/ 

)(-x<y,)= 0 

dq 

«  , 

=  2X(z/-c,-c2xi-. 

/« I 

*  *  ^(wx+l)(my+l)**'/  y i 

)(-or)-o 

^  C(mx+\){my+i) 

Rearranging  yields  the  following  set  of  ( mx  + 1  )(my  + 1)  simultaneous  equations: 


C|/?  +  Cj ^ j Xj  + . ..  +  ck  ^  \ y j  +  Xx,y,  + . . .  +  C(mx+\)(my+\)  Xj  y ,•  ^ 

ci  X xi +  c2  X xf + • .  ■ ■ + X  + ci+i  X  x/2» + •  •  • + c(„„+I)(w+i)  X  xr+vr = X  *<*/ 

C,  X  X  +  C2  X  X/X  +  •  •  ■ •  +  ck  X  T,2  +  C*+1  X  X<T,2  +  •  •  •  +  c(„u+ixm,+i)  X  *7Xyr'  =  X  Xz; 
ci  X x,x  +  c2  X  xfy- + •  •  • + ck  Xxx2  +  cA+1  X  x,2x2 + •  •  • + ^w^H^r'yr'  =  XxiTiz/ 

^Xor  +c2Xxrvr + ...+c*Xxr^r1 +ci+1ixrvr' + ... 
+w1Xm,+1)Xx^r=Xxr^ 


(48) 


The  normal  equations  represented  in  matrix  form  are: 
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Here  we  have  a  (mx  + 1  ){my  + 1)  x  (mx  +  \){my  + 1)  matrix  of  known  quantities  multiplying  an  unknown 
(mx  + 1  ){my  + 1)  x  1  coefficient  vector  to  obtain  a  known  (mx  + 1  ){my  + 1)  x  1  right-hand-side  vector.  To 
build  the  above  set  of  equations  within  a  computer  code,  one  begins  by  forming  the  n  x  (mx  + 1  ){my  + 1) 
matrix  C  and  its  (mx  + 1  ){my  + 1)  x  n  transpose  CT : 


C  = 


CT  = 


0  0 

4 l»°  ••• 

44  44 

mx  my 
x]  >\ 

0  0 
x2y2 

x\y°2  ••• 

4/2  A4l 

...  x?y? 

0  0 
x„y„ 

4y°n  ••• 

4,4n  44, 

-  OT 

x\yl 

x°2y°2 

,  ^0 

_ 1 

x\y\ 

xly°2 

•••  x„y„ 

4y\ 

4 \y\ 

44, 

44 

x\y\ 

•••  44, 

x,~yr 

mx  my 

x2  y2 

... 

n  sn  j 

and 


(50) 


Then,  the  normal  equations  defined  in  Eqs.  49  can  be  easily  developed  using  matrix  multiplication  as 
follows: 
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CTC  c  =  CTz  ,  (51) 

and  the  known  matrices  CTC  and  CTz  can  be  passed  to  a  simultaneous  equations  solver  to  solve  for  the 
unknown  coefficient  vector  c.  Note  that  when  forming  the  C  matrix  that  columns  are  formed  by  varying 
the  first  superscript  from  0  to  mx ,  then  repeating  this  pattern  incrementing  values  of  the  second 
superscript  from  0  to  m .  If  one  desires  to  omit  one  or  more  terms  from  the  fit  equation,  the 
corresponding  columns  of  the  C  matrix  and  the  corresponding  rows  of  the  CT  matrix  are  omitted  prior  to 
forming  the  CTC  and  CTz  matrices  and  the  number  of  coefficients  solved  for  is  reduced  accordingly. 

The  quantities  which  describe  the  goodness  of  the  fit  in  Eqs  8  and  9  are  for  this  case: 


-1 1/2 


a’’xy  1  n-(mx  +l)(m  +  1)£ 


EM) 


(Jy'x  1  n-{mx  +\){mv  +1)£? 


n  z 

E(z,  -  C1  -  Cj*,  -  •  •  •  -  cky,  -  ck+xx,y,  - ...  -  c(na+i){m>,+X)x",x y”,y) 


1/2 


»(52) 


and  the  correlation  coefficient,  R,  which  is  computed  from 


R  = 


~l  V2 


1  <Jl'Xy 

1  2 

CT. 


where  the  standard  deviation  of  z,  crz  = 


“1 1/2 


(53) 


Uncertainty  in  Surface  Fit  Coefficients 

The  application  of  the  uncertainty  method  requires  n  triplets  of  experimental  data  (x,  ,y, ,  z,.) ,  along  with 
associated  bias  errors,  (rJ.,  (r,,)  and  (R,). ,  and  precision  errors,  (rJ.,  (Rv).  and  (Rz)..  The  bias 
uncertainty  propagated  into  the  bfh  coefficient  ck ,  denoted  by  (Rc) ,  is  computed  from: 


wi=s( 

/=i  v 

«-!  n  ( 

+2  £  I 

f=i 

«-]  n  I 

+2Z  2 

i=l  >,‘+1 V 
n- 1  n  f 

-22  2 

i=\  j=i+ 1 V 


<?Ck) 

2 

(s)2 

dxj 

V  x/l 

dCk' 

(dck) 

d’xj 

[dXj) 

<?y,J 

<?ck] 

dck 

.dzj 

{dzj 

/=  1 


£,K 


MM, (*.),*>£  z 

7/  7  /=1  7=1+1 

MW, -22  2 


/'=!  j=i+\ 


{0ck\ 

\dz,  J 

w: 

<  a  \ 

dck 

( a  \ 
dck 

UrJ 

(<?ck ) 

( a  \ 
dck 

[dzj) 

dCk\ 

-  \ 

.dz,) 

MM 


(54) 
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where  we  have  assumed  correlated  biases  among  all  combinations  of  xf ,  y,  and  z, .  If  any  of  these 
biases  among  the  raw  data  values  can  be  shown  to  be  uncorrelated,  then  the  corresponding  terms  in 
Eq.  54  should  be  excluded. 

The  precision  uncertainty  propagated  into  the  coefficient  ck,  { ;  ,  is  obtained  from  a  similar  equation; 
however,  since  precision  errors  are  considered  to  be  random,  all  precision  errors  among  the  raw  data 
values  are  uncorrelated. 


WI-Z 

/=! 


(55) 


The  total  uncertainty,  (t/c)  ,  for  the  coefficient  ck  is  then  computed  using  either  of  the  formulas  given 
in  Eqs.  14  or  15. 
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Uncertainty  in  Predicted  Values  from  Surface  Fit 

After  calculating  the  coefficients  of  the  surface  fit  equation,  one  can  supply  new  values  X  and  Y  in  order 
to  compute  a  new  value  Z  using: 


Z(X,Y)  =  c,  +c2X  + ...  +  ctr+cMXY+ ...  + 

for  X  and  Y  chosen  such  that  xmin  <  X  <  xmax  and  >min  <Y<yn 


(56) 


The  most  general  equation  to  define  the  bias  uncertainty  that  will  propagate  into  this  fitted  Z(X,Y)  is 
found  to  be: 


A  f  dZ^ 


f  dZ^ 


/=! 
n- 1  n 


dZ; 


+  . 


w-1  n  ( 

2SI 


<?Z 


<?z 


1=1  j=i+ 1 

«-l  n 
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\dyj 

P  A  ( dZ 


f  dZ\ 

f  dz'' 

\dxj 

(  dZ\ 
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1 dy,) 

U  Zj) 
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,=1  j=i+\\V  Zj 


dZ_ 
\dz}; 


"  ( dZ ' 


(*.),(*),+*  IZV 

/•=1  j=i+\\ 


SZ: 


dZ 

\dXj, 


M,  Mj 

M,  to), 


f  dZ\ 

2 

Bl  + 

(dZ\ 

UxJ 

X 

VdY) 

B i 


+  2Z 

/=1 

n 

+  *1 


dZ\ 

f  dZ\ 

,dX) 

KdxJ 

dZ\ 

(  dZ^ 

dX) 

dzt) 

<?z) 

d Z \ 

dY) 

<?yj 

2Z 

/= 1 
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, f  az' 

(  dZ 
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Uy, 

dZ\ 

d Z \ 

U  Y) 
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(  dZ\ 

)  B,  (*,), 

(*■), 

B,  {B.\ 


(57) 


where  we  have  again  assumed  the  most  general  case  of  all  combinations  of  correlated  biases  between  X 
and  Y,  and  the  x, ,  yt  and  z, .  If  any  of  these  biases  can  be  shown  to  be  uncorrelated  then  the  appropriate 
terms  in  Eq.  57  should  be  excluded. 

The  precision  uncertainty  that  propagates  into  the  fitted  Z(X,  Y)  is  obtained  from 
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n2  JZ 

pz  =S  ~jz 


<?Z  <?Z  <?z 

To  evaluate  these  expressions,  we  need  to  determine  the  partial  derivatives:  — — ,  — —  and  along 

o  x:  o  yt  o  Zj 

3  Zd  8  Zj 

with  — —  and - .  The  derivatives  are  applied  to  Eq.  56  and  are  found  to  be: 

d  X  dY 

dZ  d  Cy  d  C-±  v  3  Ck  v  ,  3ck^  1  t  J^mxymy 


^=ZZ±  +  ^x+...  +  ^Y+^1-XY+...  +  - 
dxi  dx-,  dxy  dx,  ox ( 

m.x  my  o 

=zs  — -XpYq  9  where  k  =  q \mx  + 1)  +  p  + 1 

p= 0  <7=0  3  Xi 


8Z  8Z  ^  8Z  8Z 

with  similar  expressions  for  and  .  For  and  ^  ^  ,  we  have. 


*7  mx  my  A  7  ,  \ 

TV=SZ  PCt  Yq  and  — -  =  ££  qck  Xp  T9-1 ,  where  k  =  q  (mx +  \)  + p  +  \  (60) 

8  X  ps{  q=o  &  *  p= 0  q-\ 

Note  that  ——  and  are  functions  of  X  and  Y  and  must  be  recomputed  for  each  new  value  of 

dx-t  ’  dyt  dX 

XoxY. 

dck  dck 

Computation  of  the  expressions  in  Eqs.  54-55  and  57-59  require  the  sensitivity  derivatives:  y—  -jy 

and  — -  Using  the  method  described  earlier,  we  form  analytic  expressions  for  the  coefficients  using 

dZ; 

Cramer’s  rule  which  makes  use  of  the  CTC  and  CTz  matrices  which  are  described  in  Eqs.  49-51.  The 
expressions  for  the  coefficients  will  be  similar  in  form  to  Eqs.  37.  One  computes  derivatives  of  these 
quantities  in  a  manner  analogous  to  that  described  in  Eqs.  40-43.  The  only  change  is  in  the 
implementation  details  of  steps  3  and  7  of  the  procedure  previously  outlined;  namely,  the  replacement  of 
the  contents  of  a  column  in  a  matrix  with  its  derivative.  This  is  a  minor  change;  however,  if  one  further 
implements  the  ability  to  pick  and  choose  any  terms  in  the  desired  functional  form  of  the  surface  fit,  then 
one  must  keep  track  of  the  omitted  terms  when  incorporating  this  portion  of  the  algorithm.  In  all  other 
respects  the  calculation  of  the  sensitivity  derivatives  is  straightforward  and  identical  to  that  discussed 
previously. 

Now,  given  the  values  X  and  Y,  we  can  compute  a  value  Z  from  Eq.  56  and  an  associated  Bz  and  Pz 
from  Eqs.  57  and  58.  From  the  latter  two  quantities  we  can  determine  a  total  uncertainty  Uz  from 


22 


Eqs.  14  or  15.  This  process  should  be  repeated  for  Xj  and  Y}  which  vary  throughout  the  domain 
xmin  <  Xj  <  xmax  and  ymin<Yj<ymax  and  a  series  of  Z{X,Y)j  and  (Uz)  .  will  result.  Then,  the 
uncertainty  in  fitted  Z(X,Y)j  values  can  be  plotted  along  with  the  fitted  surface  as  follows.  Plot  the 
original  (x,  ,y,,zf)  data  values  along  with  the  fitted  surface  Z(X,Y)j  obtained  from  the  new  [Xj  ,Yj) 
data  pairs.  Then  plot  above  and  below  this  fitted  surface  the  two  additional  surfaces:  Z(X,Y)j  +{uz)j 

and  Z(X,Y)  j  -  (Uz )  .  The  latter  two  surfaces  then  give  an  indication  of  the  total  uncertainty  in  fitted 

Z(X,Y)j  values  across  the  entire  domain  xmin  <  Xj  <  xmax  and  >’min  <  Yf  <  ymax.  Since  a  global  plot  of 
this  sort  may  prove  difficult  to  extract  quantitative  information  from,  one  may  wish  to  alternatively  plot 
two-dimensional  slices  of  the  data  to  provide  local  uncertainty  behavior. 

Examples 

This  section  gives  a  few  simple  examples  of  the  uncertainty  that  one  may  expect  in  coefficients  and 
fitted  values  using  the  procedures  described  in  the  previous  sections.  The  results  provide  some  insight 
and  quantification  of  the  penalties  that  can  arise  as  a  result  of  improper  use  of  regression  techniques. 
The  approach  taken  is  to  use  a  basic  set  of  experimental  data  with  prescribed  uncertainty  and  to  perform 
a  fit  to  the  data  and  compute  the  uncertainty  propagated  through  the  fit  equations.  Then,  the  data  set  is 
altered  in  a  variety  of  ways  to  determine  the  resulting  changes  in  the  uncertainty  results.  The  data  to  be 
fitted  are  given  in  Table  1  and  include  only  precision  uncertainty  in  the  ordinate.  The  data  are  taken 
from  pp.  1 83-184  of  the  text  by  Coleman  and  Steele2  and  are  used  with  permission. 


Table  1 .  Test  data  for  examples. 


X, 

to, 

to, 

y, 

to, 

to, 

2.0 

0.0 

0.0 

2.4 

0.0 

1.0 

3.0 

0.0 

0.0 

3.0 

0.0 

1.0 

4.5 

0.0 

0.0 

3.5 

0.0 

1.0 

5.3 

0.0 

0.0 

4.5 

0.0 

1.0 

6.5 

0.0 

0.0 

4.9 

o.o 

1.0 

7.8 

0.0 

0.0 

5.6 

0.0 

1.0 

8.5 

0.0 

0.0 

6.8 

0.0 

1.0 

10.1 

0.0 

0.0 

7.3 

0.0 

1.0 

The  first  case  consists  of  a  linear  fit  to  the  data  in  Table  1.  Following  the  summary  provided  in  a 
previous  section,  the  fit  coefficients  and  the  total  uncertainties  (using  Eq.  14)  in  the  fit  coefficients  were 
computed.  The  latter  quantities  were  divided  by  the  respective  values  of  the  fit  coefficients  to  produce 
percentage  uncertainties,  and  these  data  may  be  found  in  Table  2.  For  the  given  number  of  data  points 
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and  for  the  specified  level  of  uncertainty  in  the  ordinates,  the  uncertainty  in  the  intercept  is  86%  and  in 
the  slope,  22%. 


Table  2.  Uncertainty  in  fit  coefficients. 


Type 

Cl 

C2 

c3 

Uc,/c,  (%) 

Uc2/c2  (%) 

Uc3/c3  (%) 

Case  1 

Linear 

1.0278 

0.6243 

n/a 

85.93 

21.74 

n/a 

Case  2 

Quadratic 

1.3634 

0.4857 

0.0116 

137.68 

143.51 

493.46 

Case  3 

0.6309 

n/a 

16.74 

n/a 

| 

Half  Unc 

1.0278 

0.6243 

n/a 

42.97 

10.87 

n/a 

Add  Err 

1.0278 

0.6243 

n/a 

88.20 

22.13 

n/a 

The  next  step  is  to  produce  a  series  a  fitted  values,  Y(X)jt  using  a  series  of  Xj  lying  between  the 
minimum  and  maximum  values  of  the  data  in  the  first  column  of  Table  1.  Twenty  such  values  were 
chosen,  and  since  the  original  data  values  contained  no  uncertainty  in  the  abscissa,  the  uncertainty 
associated  with  new  abscissas  was  estimated  to  be  zero.  Then,  for  each  of  these  fitted  values,  the  total 

uncertainty,  (uy).,  was  calculated  and  the  quantities  Y{X)}  +(UY).  and  Y{X)J  -(ur).  were  formed. 


Table  3.  Uncertainty  in  fitted  values. 


XN 

Yn 

UYn 

yn+uyn 

yn-uyn 

2.00 

2.28 

0.64 

2.92 

1.63 

2.43 

2.54 

0.60 

3.14 

1.95 

2.85 

2.81 

0.55 

3.36 

2.26 

1  3.28 

3.07 

0.51 

3.58 

2.57 

3.71 

3.34 

0.47 

3.81 

2.87 

4.13 

3.61  ' 

0.43 

4.04 

3.17 

4.56 

3.87 

0.40 

4.27 

3.47 

4.98 

4.14 

0.38 

4.52 

3.76 

5.41 

4.41 

0.36 

4.77 

4.04 

5.84 

4.67 

0.35 

5.03 

4.32 

6.26 

4.94 

0.36 

5.29 

4.58 

6.69 

5.20 

0.37 

5.57 

4.84 

7.12 

5.47 

0.39 

5.86 

5.08 

7.54 

5.74 

0.41 

6.15 

5.32 

7.97 

6.00 

0.45 

6.45 

5.56 

8.39 

6.27 

0.48 

6.75 

5.78 

8.82 

6.53 

0.52 

7.06 

6.01 

9.25 

6.80 

0.57 

7.37 

6.23 

9.67 

7.07 

0.62 

7.68 

6.45 

10.10 

7.33 

0.66 

8.00 

6.67 
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This  data  may  be  found  in  Table  3.  One  can  see  that  the  uncertainty  propagated  into  fitted  values  varies 
with  the  abscissa  from  a  maximum  of  0.66  near  each  end  of  the  range  to  a  minimum  of  0.35 
approximately  midway  through  the  range.  Figure  1  shows  a  plot  of  this  behavior.  Also,  contained  in 
Fig.  1  is  the  original  data  with  error  bars  representing  the  precision  uncertainty  in  the  ordinates  along 
with  the  linear  fit  satisfying  the  least  squares  criterion. 

In  addition,  quantities  which  describe  the  goodness  of  the  fit:  the  standard  error  (also  known  as  the 
standard  error  of  the  estimate,  SEE)  and  the  correlation  coefficient,  were  computed  from  Eqs.  8  and  9, 
respectively.  The  results  may  be  found  in  Table  4.  The  total  uncertainty  in  the  abscissa  for  fitted  values, 
estimated  to  be  zero,  is  listed  in  the  table  as  well  as  an  average  uncertainty  for  the  ordinate.  Finally,  for 
comparison,  the  last  column  of  Table  4  contains  the  standard  error  multiplied  by  two;  this  number  has 
often  been  used  in  the  past  as  the  uncertainty  bands  around  fitted  values  from  the  fit. 


Table  4.  Goodness  of  fit  parameters. 


Type 

R 

Std  Err 

uxN 

AvgUYN 

2  SEE 

Case  1 

Linear 

0.9863 

0.2905 

0 

0.4759 

0.5810 

Case  2 

Quadratic 

0.9848 

0.3051 

0 

0.5653 

0.6101 

Case  3 

Dble  Data 

0.9903 

0.2234 

0 

0.3574 

0.4468 

Case  4 

Half  Unc 

0.9863 

0.2905 

0 

0.2380 

0.5810 

Case  5 

Add  Err 

0.9863 

0.2905 

0.3162 

0.5303 

0.5810 

A  glance  at  Fig.  1  shows  that  the  fit  appears  to  be  reasonably  good,  and  the  fit  parameters  bear  this  out 
with  a  high  correlation  coefficient  and  a  relatively  small  standard  error  of  the  original  data  points  about 
the  fitted  line.  What  these  statistics  do  not  indicate  are  the  high  uncertainties  in  the  fit  coefficients, 
which  are  presumed  to  result  from:  the  amount  of  data  used  to  construct  the  fit  and  the  high  uncertainties 
associated  with  the  ordinates  of  the  original  data.  Fitted  values  obtained  from  this  fit  will  have  a  total 
uncertainty  associated  with  them  which  varies  from  0.66  to  0.35  or  an  average  value  of  0.48.  This  latter 
value  is  greater  than  the  standard  error,  but  less  than  the  uncertainties  in  the  ordinates  of  the  original  data 
values  and  not  as  conservative  as  the  less  accurate  2  SEE  measure  applied  in  the  past. 

Although  the  data  appear  to  be  reasonably  approximated  by  a  linear  fit,  one  might  wish  to  improve  the 
fit  by  using  a  higher  order  polynomial.  One  might  argue  that  the  coefficient  multiplying  the  quadratic 
term,  if  not  needed,  will  simply  turn  out  to  be  a  small  number,  and  the  higher  order  fit  will  do  no  harm. 
For  case  2,  a  quadratic  fit  was  computed  for  the  data  in  Table  1  using  the  same  procedure  as  described 
above.  The  coefficients  and  percent  uncertainty  for  the  coefficients  are  given  in  Table  2,  the  fit 
parameters  are  in  Table  4,  and  a  plot  of  the  results  may  be  found  in  Fig.  2.  We  see  that  for  this  particular 
data  set,  the  standard  error  and  the  correlation  coefficient  are  nearly  the  same  for  the  quadratic  and  linear 
cases.  However,  the  uncertainties  associated  with  the  coefficients,  particularly  the  quadratic  term,  are 
very  high.  The  uncertainty  in  fitted  values  varies  from  a  maximum  of  0.88  to  a  minimum  of  0.47  which 
yields  an  average  value  of  0.57,  and  this  is  higher  than  that  for  a  linear  fit  with  no  improvement  in  fit 
characteristics. 
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In  an  attempt  to  reduce  uncertainty,  the  input  data  set  listed  in  Table  1  was  approximately  doubled  in 
size,  and  a  linear  fit  was  computed  for  this  extended  data  set  for  case  3.  The  extended  data  set  was 
obtained  by  adding  an  additional  pair  of  data  values  between  each  of  the  existing  data  values  by  means 
of  linear  interpolation  in  each  coordinate,  thereby  increasing  the  data  set  from  8  to  15  data  pairs.  Each 
new  data  pair  was  assigned  zero  bias  and  precision  uncertainties  for  the  abscissas  and  a  zero  bias  and  a 
precision  uncertainty  equal  to  one  for  each  ordinate.  One  might  argue  that  the  manner  in  which  the  new 
data  w'ere  generated  may  bias  the  results,  nevertheless,  the  improvement  was  surprisingly  meager.  The 
coefficients  and  percent  uncertainty  for  the  coefficients  are  given  in  Table  2,  and  the  fit  parameters  are  in 
Table  4.  The  computed  coefficients  varied  slightly  due  to  the  change  in  the  input  data,  but  were  close 
enough  to  case  1  to  afford  a  reasonable  comparison.  There  was  some  improvement  in  the  uncertainty  in 
the  fit  coefficients:  a  reduction  from  86%  to  69%  for  the  intercept  and  a  reduction  from  22%  to  17%  for 
the  slope.  Similarly,  the  uncertainty  in  fitted  values  varied  from  a  maximum  of  0.51  to  a  minimum  of 
0.26  for  an  average  value  of  0.36.  The  correlation  coefficient  improved  from  0.986  to  0.99  and  the 
standard  error  decreased  from  0.29  to  0.22. 

The  fourth  example  was  designed  to  test  the  effect  of  halving  the  uncertainty  in  the  input  data.  Thus,  the 
eight  data  pairs  listed  in  Table  1  were  used  with  the  exception  that  the  precision  uncertainties  for  the 
ordinates  were  reduced  to  0.5.  The  results  of  the  linear  fit  are  listed  in  the  appropriate  tables,  and  we 
find  a  marked  decrease  for  this  case  when  compared  to  case  1 .  The  uncertainty  in  the  intercept  was 
reduced  from  86%  to  43%  and  for  the  slope  from  22%  to  11%.  The  uncertainty  in  fitted  values  varied 
from  a  maximum  of  0.33  to  a  minimum  of  0.18  for  an  average  value  of  0.24.  The  correlation  coefficient 
and  standard  error  were  identical  to  that  of  the  first  case  because  these  parameters  are  formed  from  the 
input  coordinate  pairs  which  did  not  change.  Notice  the  considerable  difference  between  the  average 
uncertainty  in  the  ordinate,  0.24,  and  the  2  SEE  parameter  which  remained  the  same  as  for  case  1  at 
0.58;  the  2  SEE  parameter  is  an  excessively  conservative  bound  for  this  case.  For  a  given  fit  order,  a 
reduction  in  uncertainty  in  the  data  to  be  fitted  is  far  more  effective  in  reducing  uncertainty  associated 
with  fit  coefficients  and  with  fitted  values  than  simply  increasing  the  amount  of  data  to  be  fitted. 


Table  5.  Test  data  for  case  5. 


*/ 

M, 

00, 

y, 

W, 

M, 

2.0 

0.1 

0.3 

2.4 

0.1 

1.0 

3.0 

0.1 

0.3 

3.0 

0.1 

1.0 

4.5 

0.1 

0.3 

3.5 

0.1 

1.0 

5.3 

0.1 

0.3 

4.5 

0.1 

1.0 

6.5 

0.1 

0.3 

4.9 

0.1 

1.0 

7.8 

0.1 

5.6 

0.1 

1.0 

8.5 

0.1 

0.3 

6.8 

0.1 

1.0 

10.1 

0.1 

0.3 

7.3 

0.1 

1.0 

Finally,  for  the  fifth  case,  reasonable  bias  and  precision  errors  were  chosen  for  both  the  abscissas  and  the 
ordinates  for  the  data  listed  in  Table  1  in  order  to  better  simulate  the  case  of  experimentally  derived  data. 
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Although,  in  general,  these  errors  may  vary  from  data  pair  to  data  pair  in  the  set,  for  simplicity,  the  same 
estimates  were  used  for  each  data  pair.  The  data  set  used  for  this  example  may  be  found  in  Table  5. 

In  order  to  calculate  the  bias  error  propagated  into  the  coefficients  (Bc)k  from  Eq.  12,  one  must  consider 
which  biases,  if  any,  are  correlated.  Since  the  same  numerical  values  were  prescribed  for  each  of  the 
abscissas  and  ordinates,  a  reasonable  assumption  is  that  the  (#,).  are  correlated  with  each  other  and  the 

[By )  are  correlated  with  each  other,  but  that  the  (Bx).  are  not  correlated  with  the  (£y)  .  For  calculating 

the  bias  error  propagated  into  fitted  values  By  from  Eq.  31,  one  must  estimate  the  expected  bias  in  the 
abscissa  Bx  for  a  fitted  data  pair  and  decide  if  this  bias  is  correlated  with  any  other  biases.  Bx  was 

assigned  values  of  0.1  for  each  new  data  pair  and  was  assumed  to  be  correlated  with  the  (Bx).  but  not 

with  the  (i^).  The  results  are  listed  in  the  appropriate  tables.  The  computed  fit  coefficients  are 

identical  to  case  1  since  the  (x,  ,y(  )  data  pairs  are  unchanged,  and  the  uncertainty  propagated  into  the  fit 
coefficients  is  about  the  same.  The  correlation  coefficient  and  standard  error  are  unchanged  from  case  1 , 
and  the  total  uncertainty  in  the  abscissa  for  fitted  values  is  0.32  which  is  computed  using  Eq.  14.  The 
uncertainty  in  fitted  ordinates  varied  from  a  maximum  of  0.71  to  a  minimum  of  0.42  for  an  average 
value  of  0.53  which  is  about  10%  higher  than  for  case  1.  The  2  SEE  bound  remained  unchanged.  The 
fact  that  this  case  differs  little  from  case  1  indicates  that  the  uncertainty  calculations  are  dominated  by 
the  precision  uncertainty  in  the  ordinates  of  the  input  data  and  that  the  largest  improvement  will  result 
from  a  decrease  in  this  uncertainty  as  shown  by  case  4. 
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Conclusions 

This  paper  has  reviewed  the  manner  in  which  curve  fits  and  surface  fits  of  arbitrary  order  which  satisfy 
the  least  squares  criterion  are  constructed.  The  method  for  the  computation  of  uncertainty  propagated 
into  fit  coefficients  and  into  fitted  values  obtained  from  the  fit  for  both  curve  and  surface  fits  has  been 
outlined.  The  primary  obstacle  to  the  efficient  implementation  of  these  calculations:  the  determination 
of  the  sensitivity  derivatives,  has  been  removed  with  the  explanation  of  an  analytic  method  for  the 
calculation  of  these  derivatives  for  arbitrary  order  curve  and  surface  fits.  A  detailed  prescription 
describing  one  possible  approach  for  the  implementation  of  the  method  was  provided.  To  illustrate  the 
power  of  the  techniques,  simple  examples  were  chosen  and  contrasted  using  a  generic  data  set  which  is 
typical  of  experimentally  derived  data.  The  results  showed  the  general  behavior  of  the  calculations 
under  a  variety  of  conditions.  In  particular,  the  danger  of  fitting  a  higher  order  polynomial  to  a  data  set 
when  such  a  model  is  not  warranted,  was  illustrated.  For  this  particular  data  set,  a  reduction  in  the 
uncertainty  associated  with  the  input  data  was  more  effective  at  producing  reductions  in  uncertainty 
associated  with  fit  coefficients  and  fitted  values  than  simply  using  more  data  with  the  original 
uncertainty  levels.  Although  the  equations  are  somewhat  tedious,  these  methods  may  be  readily 
implemented  within  a  computer  program  making  the  computation  of  uncertainty  for  fits  employing  the 
method  of  Least  Squares  a  routine  matter. 
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Appendix 
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Testing  the  Surface  Fit  Implementation 

To  test  a  surface  fitting  computer  code,  one  commonly  generates  a  set  of  data  triplets,  (x,  ,yjf  z;) ,  which 
satisfy  a  known  function  z-  f  (x,y)  formed  from  the  product  of  two  polynomials.  Then,  the  data  are 
input  into  the  program,  the  user  selects  mx  and  my  to  match  the  characteristics  of  the  known  function, 
and  the  output  will  be  a  set  of  coefficients  which  should  match  those  of  the  function  which  generated  the 
data.  Specifically,  a  set  of  (x,.  ,y,)  are  produced  in  some  convenient  manner  and  then  input  to  the 
function  z  =  f(x,y)  to  obtain  z,  and  complete  the  data  triplets  (x,.  .  This  seemingly  innocuous 

procedure  can  lead  to  frustration  because  the  generation  of  test  data,  (x,  ,y,) ,  for  input  into  z  =  /  (x,y) 
cannot  be  chosen  arbitrarily  if  one  desires  a  unique  solution!  To  see  how  problems  can  arise,  consider 
the  normal  equations,  Eqs.  49,  for  the  case  mx  =  my  =  1  which  become: 


I*,1*0 
Xx't,°  X*/V,° 
X^1 
.X*,y 


X*W  £xiy! 

Xx!y!  ’YjXfy] 
£xly? 
Zx'y,2  Yjxfyf 


V 

x*,% 

x*;a 

Xx,% 

-C4_ 

M 

vT 

jnT* 

i _ 

(Al) 


Now  if ,  for  convenience,  one  were  to  choose  x,  =  y,  yielding  data  triplets  of  the  form  (x/  ,x,.,z/),  then 
when  developing  the  coefficient  matrix  above,  row  2  and  row  3  would  become  identical: 


X*/  X*,2  X*,2  X*-3  -row  2 
£x#*  Zxf  Xx<3  -row  3 

This  causes  the  matrix  to  become  singular  with  zero  determinant.  In  other  words,  we  have  three 
independent  equations  in  the  four  unknowns  yielding  a  one-parameter  family  of  solutions  instead  of  a 
single  unique  solution.  More  generally,  if  one  chooses  y,  =  ax,.  +  b ,  one  can  show  that 
row  3  =  tf(row  2)  +  b( row  l)  which  will  again  lead  to  a  non-unique  solution  to  the  system  of  equations. 
The  point  is  that  the  method  will  fail  to  produce  a  unique  solution  when  the  x,  and  y,  are  chosen  such 
that  they  are  linearly  dependent.  The  data  to  be  fitted  must  satisfy  the  requirement  that  z  =  / (x,_y)  be  a 
function  of  two  linearly  independent  variables.  This  situation  is  unlikely  to  occur  when  fitting 
experimentally  derived  data  and  is  usually  only  encountered  when  generating  artificial  data  for  code 
validation. 
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Table  A1 .  Test  data  generation  scheme. 
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1 

3 

2 

1 
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2 

1 

2 

2 

5 
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2 

2 

3 
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2 

3 

3 

1 

7 

3 

1 

3 

2 

8 

3 

2 

3 

3 

9 

3 

3 

One  simple  scheme  for  creating  test  data  is  to  choose  an  integer,  /,  such  that  the  desired  number  of  data 
triplets  will  be  n  =  l2;  then  using  two  loops,  with  indices  r  and  s,  one  constructs  the  x,  and  y,  using: 

x;  =rx0  and  y,  =syQ  where  i  =  {r-\)l  +  s 
1  <r<l,  1  <s<l,  1  <i<n;  x0  and  y0  are  arbitrary  scale  factors  ' 

For  example,  if  one  chooses  n  =  9  and  xQ=y0=  1.0,  this  scheme  produces  the  data  shown  in  Table  1 . 
One  then  chooses  an  appropriate  z  =  f(x,y )  of  the  form 


mx  my 

xpyq ,  where  k  =  q  ( mx  + 1)  +  p  + 1  and  1  <  ^  <  (mx  + 1  )(my  + 1)  ,  (A4) 

p- 0  9=0 

and  specifies  mx  and  mv .  Finally,  the  coefficients  for  this  function  must  be  chosen.  To  avoid  large 
numbers  in  the  coefficient  matrix  for  increasingly  large  choices  of  mx  and  my ,  one  can  define  the 
coefficients  to  be 


c-ir 

k 


(A5) 


Summarizing,  one  generates  n  pairs  of  x,  and  y,  using  Eq.  A3,  then  inserts  these  values  in  Eq.  A4  to 
obtain  the  zt  to  complete  each  data  triplet.  When  this  test  data  is  input  into  the  surface  fitting  code,  it 
should  return  the  fit  coefficients  defined  in  Eq.  A5  with  appropriate  choices  for  mx  and  my . 
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